Deconstructing classical water models at interfaces and in bulk 
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Using concepts from perturbation and local molecular field theories of liquids we divide the po- 
tential of the SPC/E water model into short and long ranged parts. The short ranged parts define 
a minimal reference network model that captures very well the structure of the local hydrogen 
bond network in bulk water while ignoring effects of the remaining long ranged interactions. This 
deconstruction can provide insight into the different roles that the local hydrogen bond network, 
dispersion forces, and long ranged dipolar interactions play in determining a variety of properties 
of SPC/E and related classical models of water. Here we focus on the anomalous behavior of the 
internal pressure and the temperature dependence of the density of bulk water. We further utilize 
these short ranged models along with local molecular field theory to quantify the infiuence of these 
interactions on the structure of hydrophobic interfaces and the crossover from small to large scale 
hydration behavior. The implications of our findings for theories of hydrophobicity and possible 
refinements of classical water models are also discussed. 



I. INTRODUCTION 

Classical empirical water potentials involving fixed 
point charges and Lennard- Jones (LJ) interactions were 
introduced in the first computer simulations of water 
forty years ago and modern versions are widely used even 
today in many biomolecular and materials-based simula- 
tions. Two recent reviews [ij, [331 have focused on this 
wide class of model potentials and assessed their per- 
formance for a broad range of different structural and 
thermodynamic properties, some of which were used as 
targets in the initial parameterization of the models. De- 
spite known limitations associated with the lack of molec- 
ular flexibility and polarizability, they qualitatively and 
often quantitatively capture a large number of properties 
of water and often represent a useful compromise between 
physical realism and computational tractability. 

Given the simple functional forms of the intermolecular 
potentials it may seem surprising that such good agree- 
ment is possible. But recent work has shown that even 
simpler models where particles interact via isotropic re- 
pulsive potentials with two distinct length scales are able 
to qualitatively reproduce certain characteristic dynamic 
and thermodynamic anomalies of bulk water [8|, |42| . |43| . 
Similarly in dense uniform simple liquids a hard-sphere- 
like repulsive force reference system can give a good de- 
scription of the liquid structure, and this in turn permits 
thermodynamic properties to be determined by a simple 
perturbation theory [1^ \^ . 

This suggests it should be useful to analyze the con- 
struction and predictions of empirical water potentials 
from the perspective of perturbation theory of unifornr 
fluids and the related Local Molecular Field (LMF) the- 
ory [H, [H [H, HO, [M EO] • LMF theory provides a more 
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general approach applicable to both uniform and nonuni- 
form fluids and gives strong support to the basic idea of 
perturbation theory that in a uniform fluid slowly vary- 
ing long ranged parts of the intermolecular interactions 
have little effect on the local liquid structure. 

To apply these ideas to water we divide the intermolec- 
ular interactions in a given water model into appropri- 
ately chosen short and long ranged parts. In this con- 
text, it is conceptually useful to consider separately the 
slowly varying long ranged parts of both the LJ inter- 
actions and the Coulonrb interactions. This deconstruc- 
tion of the water potential via LMF theory provides a 
hierarchical framework for assessing separately the con- 
tributions of (i) strong short ranged interactions leading 
to the local hydrogen bonding network, (ii) dispersive at- 
tractions between water molecules, and (iii) long ranged 
dipolar interactions between molecules. Disentangling 
these contributions without the insight of LMF theory 
is very difficult due to the multiple contributions of the 
point charges and the LJ interactions in standard molec- 
ular water models 

In uniform systems, the long ranged forces on a given 
water molecule from more distant neighbors tend to can- 
cel 40, 41). The remaining strong short ranged forces be- 
tween nearest neighbors arise from the interplay between 
the repulsive LJ core forces and the short ranged attrac- 
tive Coulomb forces between donor and acceptor charges. 
These forces determine a minimal reference model that 
can accurately describe the local liquid structure - the 
hydrogen-bond network for bulk water. The slowly vary- 
ing parts of the intermolecular interactions are not im- 
portant for this local structure and could be varied es- 
sentially independently to help in the determination of 
other properties as is implicitly done in the full irrodel. 
Based on previous work with LMF theory [28|, 29, 3l|, 
we examine two basic areas where we expect the differ- 
ent contributions to play varying but important roles - 
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FIG. 2. (Color online) Optimal hydrogen bonding config- 
uration of water taken from two molecules in ice Ih. LJ 
cores are depicted as gray transparent spheres with a diam- 
eter <tlj = 3.16 A, while the hydrogen bond between waters 
with oxygens separated by 2.75 A is illustrated by a dashed, 
blue cylinder. Oxygen and hydrogen atoms are colored red 
and white, respectively. 



FIG. 1. (Color online) Schematic diagram of the SPC/E water 
model listing its various geometric parameters and interaction 
parameters. The O-H bond length and H-O-H angle are fixed, 
such that the molecule is rigid. The LJ well depth is eLj ~ 
0.65 kj/mol. The oxygen site is depicted as a large red circle, 
while the hydrogen atoms are shown as smaller, gray circles. 



bulk thermodynamics and nonuniform structure. The 
short ranged interactions responsible for the hydrogen- 
bonding network are clearly necessary in all cases. LMF 
theory allows us to determine the relative importance 
of dispersive attractions and long-ranged dipolar attrac- 
tions in these applications using simple analytical correc- 
tions for thermodynamics and an cfFcctivc external field 
for nonuniform structure. 

In the next section wc discuss the separation of the wa- 
ter potential into short and long ranged parts and show 
that a minimal short ranged reference model can very 
accurately describe atom-atom correlation functions and 
other properties of the hydrogen bond network in bulk 
water. In Sec. 3 we examine the thermodynamic impli- 
cations of this picture, focusing on two anomalous prop- 
erties of bulk water qualitatively well described by the 
full water model, the density maximum at one atmo- 
sphere pressure and the behavior of the "internal pres- 
sure" {dE/dV)T, which has a temperature and density 
dependence opposite to that of most organic solvents. 
Then in Sec. 4 we look at the effects of the unbalanced 
long-ranged electrostatic and dispersion forces at aque- 
ous interfaces. LMF theory reveals strong similarities 
between the behavior of water at the liquid-vapor in- 
terface and a planar hydrophobic wall, consistent with 
previous work, and provides new insight into the relative 
importance of electrostatic and dispersion forces and the 
transition from small to large scale hydrophobicity as the 
radius of a nonpolar solute is increased. 



II. LOCAL HYDROGEN BONDS IN FULL AND 
TRUNCATED WATER POTENTIALS 



In this paper we consider one of the simplest and most 
widely used water models, the extended simple point 
charge (SPC/E) model Q, but similar ideas and conclu- 
sions apply immediately to most other members of this 
class. As shown in Fig. [I] SPC/E water consists of a LJ 
potential as well as a negative point charge centered at 
the oxygen site. Positive point charges are fixed at hy- 
drogen sites displaced from the center at a distance of 1 A 
with a tetrahedral HOH bond angle. It is a remarkable 
fact that this simple model can reproduce many struc- 
tural, thermodynamic, and dielectric properties of bulk 
water as well as those of water in nonuniform environ- 
ments around a variety of solutes and at the liquid- vapor 
interface. 

In the following we use the perspective of perturbation 
and LMF theory to help us see how this comes about. We 
use these ideas here not to suggest more efficient simu- 
lations using short ranged model potentials but rather 
as a method of analysis that provides physical insight 
into features of the full model as well. Since a detailed 
description and justification of LMF theory is given else- 
where [30|, we will focus on qualitative arguments and 
just quote specific results when needed. 

Fig. [5] gives some insight into why a perturbation pic- 
ture based on the dominance of strong short ranged forces 
in uniform environments could be especially accurate for 
bulk SPC/E and related water models. This shows two 
adjacent water molecules with a separation of 2.75 A that 
form an optimal hydrogen bond as seen in the structure 
of ice Ih. Hydrogen bonding in this model is driven by 
the very strong attractive force between opposite charges 
on the hydrogen and oxygen sites of adjacent properly 
oriented molecules. Proper orientation permits similar 



strong bonds to form with other molecules, leading to 
a tetrahedral network in bulk water. The gray circles 
drawn to scale depict the repulsive LJ core size as de- 
scribed by the usual parameter aLj = 3.16 A. The sub- 
stantial overlap indicates a large repulsive core force op- 
posing the strong electrostatic attraction, finally result- 
ing in a nearest neighbor maximum in the the equilibrium 
oxygen-oxygen correlation function of 2.75 A. 

It is interesting to note that the first BNS water 
model introduced in 1970 used a smaller core size glj — 
2.82 A [J]. However a much larger LJ core with strong 
core overlap at typical hydrogen-bond distances is a 
common property of almost every water model intro- 
duced since then and seems to be a key feature needed 
to get generally accurate results from simple classical 
point charge models. Evidentially the highly fluctuat- 
ing local hydrogen-bond network in these models arises 
from geometrically-frustrated "charge pairing" , where 
the strong LJ core repulsions and the presence of other 
charges on the acceptor water molecule oppose the close 
approach of the strongly-coupled donor and acceptor 
charges. 

We can test the accuracy of this picture by considering 
various truncated or "short" water models where slowly 
varying long ranged parts of the Coulomb and LJ inter- 
actions in SPC/E water arc completely neglected. We 
first consider a Gaussian-truncated (GT) water model, 
aheady studied by LMF theory [21 [gl [3l|. Here 
the Coulomb potential is separated into short and long 
ranged parts as 

, , 1 erfc(r/cr) erf(r/cr) , , , , ,_,, 

r r r 

where erf and crfc are the error function and comple- 
mentary error function, respectively. The short-ranged 
VQ{r) is the screened electrostatic potential resulting from 
a point charge surrounded by a neutralizing Gaussian 
charge distribution of width a. Hence wo(?') vanishes 
rapidly at distances r much greater than a while at dis- 
tances less than a the force from VQ{r) approaches that 
of the full 1/r potential. 

In GT water, depicted in Fig. [3^, the Coulomb poten- 
tial associated with each charged site in SPC/E water is 
replaced by the short-ranged vq with no change in the LJ 
interaction. As suggested by Fig. [U important features 
of the local hydrogen-bond network should be well cap- 
tured by such a truncated model if the cutoff distance 
controlled by the length parameter a in Eq. ([T]) is cho- 
sen larger than the hydrogen bond distance. Following 
Rcfs. 2^ and [31|, here we make a relatively conservative 
choice of cr = 4.5 A, but values as small as 3 A give es- 
sentially the same results. The circles are drawn to scale 
with diameters a and ctlj. 

The basic competition between very strong short 
ranged repulsive and attractive forces in the hydrogen 
bond depicted in Fig. [5] should be captured nearly as 
well by an even simpler reference model where the LJ 
potential is truncated as well, and replaced by the re- 
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FIG. 3. (Color online) Schematic diagrams of (a) GT and (b) 
GTRC water models. Truncated interactions are indicated 
by dashed lines, while full interaction potentials are indicated 
by solid lines. LJ interactions are represented by black lines, 
while oxygen and hydrogen electrostatic interaction potentials 
are shown as red and gray lines, respectively. 



pulsive force reference potential uo(^) used in the WCA 
perturbation theory for the LJ fluid [401 ■ The result- 
ing Gaussian truncated repulsive core (GTRC) model is 
schematically shown in Fig. [5]d. 

As discussed in perturbation theories of simple liq- 
uids [la, [40] , a well-chosen reference system should accu- 
rately reproduce bulk structure present in the full system 
at the same fixed density and temperature. As illustrated 
by the pair distribution functions in Fig. 01 bulk GT and 
GTRC water models have a liquid state structure virtu- 
ally identical to that in the full SPC/E model. This very 
good agreement is also reflected in other properties of 
the hydrogen-bond network. We directly examined the 
hydrogen bonding capabilities of GT and GTRC water 
models through the calculation of the average number of 
hydrogen bonds per water molecules, {uhb), as well as 
the probability distribution of a water molecule taking 
part in uhb hydrogen bonds, P (nHs), using a standard 
distance criterion of hydrogen bonds, Rqo < 3.5 A and 
feoo' < 30°, where Rqo is the oxygen-oxygen distance 
and 6'hoo' is the angle formed by the H-0 bond vector 
on the hydrogen bond donating water molecule and the 
0-0' vector between the oxygen on the donor water (O) 
and the acceptor oxygen (O') (23. For both GT and 
GTRC water models, {uhb) and P^hhb) were calcu- 
lated at temperatures ranging from 220-300 K, and were 
found to be nearly identical to the analogous quantities 
in the full SPC/E model (data not shown). These find- 
ings give credence to the idea that these two truncated 
models reproduce the hydrogen-bond network of the full 
model to a high degree of accuracy. 

These truncated models offer a minimal structural rep- 
resentation of bulk water as a fluctuating network of 
short ranged bonds determined mainly by the balance 
between the very strong electrostatic attraction between 
donor and acceptor charges and the very strong repul- 
sion of the overlapping LJ cores. We can view them as 
primitive water models in their own right, analogous to 
other simplified models recently proposed, which capture 
very well arguably the most important structural fea- 




FIG. 4. (Color online) (a) Oxygen-oxygen, oxygen-hydrogen, and hydrogen- hydrogen site-site pair distribution functions, 
gooir), goH{r), and gHH{r), respectively, for the three water models under study at T = 300 K and v = 30.148 A'', gnn and 
goH have been shifted by 0.5 and 1 units, respectively, for clarity, (b) Differences between gooir) of the full model and that 
of the designated reference systems, Agoo{r)- 



ture of bulk water, the hydrogen bond network, and it 
is instructive to see what other properties of water such 
minimal network models can describe. But corrections 
from neglected parts of the intermolecular interactions 
are certainly needed for bulk thermodynamic and dielec- 
tric properties and for both structure and thermodynam- 
ics of water in nonuniform environments. LMF theory 
provides a more general framework where the truncated 
models are viewed as useful reference systems that can 
be systematically corrected to achieve good agreement 
with full water models. We will use both viewpoints in 
the next section. 



III. THERMODYNAMIC ANOMALIES OF 

BULK WATER USING FULL AND TRUNCATED 

WATER POTENTIALS 

A. Density mEtximum 

Now we turn our attention to the thermodynamics of 
bulk water. For a fixed volume V, temperature T, and 
number of molecules N, the pressure and other thermo- 
dynamic properties of the GT and GTRC systems will 
not generally equal those of the full system. However, 
because of the accurate reference structure, we can cor- 
rect the thermodynamics using simple mean-field (MP) 
arguments. Thus we can define the pressure in the full 
system to be the sum of the short-ranged reference pres- 
sure and a long-ranged correction, P = Pq + Pi. 

Simple corrections to the energy and pressure of 
the GT model from this perspective were recently de- 
rived [31|. With a = 4.5 A, these corrections are rela- 
tively small and were ignored in most earlier work using 
truncated water models but they are conceptually impor- 



tant in revealing the connections between truncated mod- 
els and perturbation theory and are required for quanti- 
tative agreement. The pressure correction Pi ~ Pf for 
the GT model arises only from long-ranged Coulomb in- 
teractions and is given as 
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where e is the dielectric constant. 

In the case of the GTRC model, the need for a thermo- 
dynamic correction is much more obvious since we have 
to correct for the absence of LJ attractions as well. Here 
we adopt the simple analytic correction used in the van 
der Waals (vdW) equation derived from WCA theory for 



the LJ fluid, as discussed in Ref. [3£ 
for the GTRC potential, where 
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and ui is attractive part of the LJ potential. This sim- 
ple approximation does not give quantitative results but 
does capture the main qualitative features and we use it 
here to emphasize the point that both the long ranged 
Coulomb and dispersion force corrections to bulk GTRC 
water can be treated by simple perturbation methods. 

We can test the accuracy of these corrections by using 
them to help determine the temperature Tmd at which 
the density maximum of the full SPC/E water model at 
a constant pressure of 1 atm should occur. This can 
alternatively be defined as the temperature at which the 
thermal expansion coefficient, ap, is zero. Accordingly, 
we seek to evaluate ap using the relation 

^P^-(^] =--(S) (Ih) ' (4) 
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FIG. 5. (Color online) The dependence of (a) density and (b) 
pressure for the corrected reference models as a function of 
temperature. The analogous quantities for the models with- 
out MF corrections are shown in (c) and (d), respectively. 
p{T) is calculated at a constant pressure of 1.0 atm and P{T) 
is calculated at a fixed volume of v = 30.148 A^. Full SPC/E 
data for p{T) at constant P was taken from the work of Ash- 
baugh et al. [3|. 



where v = V/N is the volume per particle. Using 
the last expression we can determine where the quan- 
tity (dP/dT),^ = by monitoring the corrected pres- 
sure of the reference systems while changing the tem- 
perature at a fixed density. This can be done by simula- 
tion in the canonical ensemble. The fixed density ensures 
that the structure of the reference and full systems are 
very similar, as assumed in the derivation of the correc- 
tions in Eqs. ([2]) and ([3]). We can also determine Tmd 
through the first equality in Eq. (|4]) by finding where 
{dv/dT)p = 0. Thus we simulate GT and GTRC water 
at constant pressures of Pq = P — Pi, where P = 1.0 atm 
is the pressure in the full system. Note that the correc- 
tion P' = P^{T;e{T)) is temperature-dependent, as is 
the dielectric constant e, so that we are not moving along 
an isobar in Pq, but an isobar in P. The temperature- 
dependent values of e were taken to be the experimental 
values (iq . 

Figs. [5^ and [5)d give the density and pressure of full 
SPG/E water and the corrected reference models as a 
function of temperature. As expected, the inclusion of P' 
in the pressure of GT water quantitatively corrects the 
density and pressure of this system. However, the MF 
correction applied to GTRC water, Pq = P — Pj' -I- ap^, 
is not as accurate, although the dependence of p on T is 
qualitatively well captured. These remaining errors arise 
from our use of the simple van der Waals ap"^ correction 
for the long ranged part of the LJ potential. This level of 
agreement is typical when this correction is used in pure 
LJ fluids [33 and a full WCA perturbative treatment of 
the attractive portion of the LJ potential in GTRC water 
would likely lead to quantitatively accurate results [401 • 

We now turn to the alternate and less accurate inter- 



pretation of the GT and GTRC models as primitive water 
models in their own right. Do these models at an uncor- 
rected pressure of 1 atm have a density maximum and 
how well does it compare to that of the full model? To 
that end, we find where {dv/dT)p = in each model by 
varying the temperature along an isobar using MD simu- 
lations in the isothermal-isobaric ensemble at a constant 
pressure of 1 atm. By requiring the same pressure in the 
full and reference models, we probe structurally differ- 
ent state points in general and there is no guarantee that 
the density and temperature of the reference systems at 
a density maxima (if present) will be similar to that in 
the full system. Nevertheless Fig. [5]; shows that the GT 
model does have a density maximum very similar to that 
of the full model. This is because the pressure correction 
to the density from the long-ranged Coulomb interactions 
in Eq. ([2|) is very small on the scale of the graph. In con- 
trast, the uncorrected GTRC model does not exhibit a 
density maximum at P = 1.0 atm, even upon cooling to 
50 K (data not shown). 

These results should be compared to earlier work where 
the TIP4P water potential was approximated by a sim- 
pler "primitive model" [20[. In that work, the repul- 
sive LJ core was mapped onto a hard-sphere potential, 
hydrogen bonding was captured by a square-well poten- 
tial, and long-ranged dipole-dipole interactions were rep- 
resented with a dipolar potential. The equation of state 
was found using a perturbative approach, and thermody- 
namic quantities were analyzed. The authors of Ref. [20| 
found that the inclusion of dispersion forces does not 
lead to a density maximum, and only when both dis- 
persive interactions and long-ranged dipole-dipole inter- 
actions were taken into account did a density maximum 
appear. 

To provide some understanding of these differing re- 
sults, we analyze the structure of the uncorrected GT 
and GTRC reference models in comparison to the full 
model at the common pressure of one atmosphere. The 
oxygen-oxygen radial distribution functions, 500(^)1 for 
each of the three water models at P = 300 K are de- 
picted in Fig. [6^. The GT model is in good agreement 
with the full model, consistent with its accurate descrip- 
tion of the bulk water density and the density maximum. 
In contrast, as shown later in Fig.|9l the coexisting liquid 
density of GTRC water is about about 15% lower than 
that of the full water model. Nevertheless the first peak 
of 500 (^) in GTRC water is higher than that of the full 
water model due to better formation of local hydrogen 
bonds. As shown in the inset, a molecule of GTRC wa- 
ter has slightly fewer hydrogen bonds on average than 
full and GT water models for temperatures higher than 
240 K. However the hydrogen bond efficiency shown in 
Fig. [6b, 
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where (tinn) is the average number of nearest-neighbors 
satisfying Poo < 3.5 A, indicates that GTRC water is 
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FIG. 6. (Color online) (a) The oxygen-oxygen pair distribution function, goo{r), for the three water models at T = 300 
K. Inset: The average number of hydrogen bonds per water molecule as a function of temperature, (nHB(T)), for full and 
truncated water models, (b) Hydrogen bonding efficiency t^hb as a function of temperature. All results were obtained at a 
constant uncorrected pressure of 1 atm. 



about 10 percent more efficiently hydrogen bonded to its 
available neighbors at aU temperatures. In this sense the 
low density GTRC water at P = 1.0 atm is structurally 
more icc-likc than the full water model. 

These results provide some insight into earlier first 
principles simulations o f liq uid water using density func- 
tional theory [2l|, [S^, [S^- The standard exchange- 
correlation functionals used there can give a good de- 
scription of local hydrogen bonding, but do not include 
effects of van der Waals interactions. These simulations 
produced a decrease in the bulk density of water accom- 
panied by increased local structural order very similar 
to that seen here for GTRC water. Moreover, when 
dispersive interactions were crudely accounted for, they 
observed much better agreement with experiment, in 
complete agreement with our findings for perturbation- 
corrected GTRC water. 

Our results indicate that van der Waals attractions 
play the role of a cohesive energy needed to achieve the 
high density present in SPC/E water at low pressure, as 
demonstrated by the qualitative accuracy of Eq. ([3]) and 
the good agreement of the GT model. Evidentially a den- 
sity maximum can arise only when additional somewhat 
less favorably bonded molecules are incorporated into the 
GTRC network to produce the full water density. If the 
local hydrogen bond network of water at the correct bulk 
density is properly described, long ranged dipolar forces 
arc not needed to obtain the correct behavior of p{T). 
Indeed, LJ attractions are not needed either provided 
that the proper bulk density is prescribed by some other 
means. Thus we found that if GTRC water is kept at a 
high constant pressure of 3 katm, where its bulk density 
is close to that of the full water model at ambient con- 
ditions, a density maximum is also observed (data not 
shown) . 



B. Internal pressure 

We further employ the reference water models to ex- 
plain the anomalous "internal pressure" of water [33|. 
For a typical van der Waals liquid, the internal pressure 
is given by Pi = {d£/dv)rp « ap^ for low to moderate 
densities, where e = E/N is the energy per molecule. In 
fact, it was recently shown by computer simulation that 
the portion of the internal pressure due to the attractions 
in a LJ fluid displays this ap^ dependence even at high 
densities [1^. Water, on the other hand, displays a neg- 
ative dependence of Pi on density. It is this anomalous 
behavior that we seek to explain. 

We begin by partitioning the internal energy of the 
system as 

e = e^-^ + e«, (6) 

where e'""' is the Lennard- Jones contribution to the en- 
ergy and £* is the energy due to charge-charge interac- 
tions (note that the change in kinetic energy when per- 
turbing the volume at constant T is zero, so we only 
consider the potential energy). We can then write the 
internal pressure as 
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This decomposition of Pi will allow us to determine which 
molecular interactions arc responsible for the strange de- 
pendence of this quantity on p. 

Fig. [2] suggests the following qualitative picture. At 
a given temperature and density the dominant hydrogen 
bond contribution to the energy e is determined from 
the balance between strong repulsive forces from the LJ 
cores and strong attractions from the more slowly vary- 
ing Coulomb interactions between donor and acceptor 
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FIG. 7. (Color online) (a) The electrostatic contribution to 
the internal pressure, P^, and (b) the analogous contribution 
from LJ interactions, Pi . The total internal pressure as a 
function of density is shown in the inset. Lines are guides to 
the eye. 
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FIG. 8. (Color online) (a) The electrostatic contribution to 
the internal pressure, Pf , and (b) the analogous contribution 
from LJ interactions, Pt^ , both as a function of temperature 
at a fixed volume of i; = 29.9 A'^. The total internal pressure 
as a function of temperature is shown in the inset. Lines are 
guides to the eye. 



charges. The Coulomb contribution P^ to the internal 
pressure Pi (T, p) is positive since a small positive change 
in volume reduces the negative Coulomb energy and sim- 
ilarly the LJ core contribution to P^"^ is negative. If the 
density is now varied at constant temperature we would 
expect the changes in Pi{T,p) to be dominated by the 
rapidly varying LJ core forces. 

Conversely, to the extent that the repulsive LJ cores 
are like hard spheres, they would contribute no tempera- 
ture dependence to the internal pressure at fixed density. 
Thus we expect the more slowly varying Coulomb forces 
to largely determine how the internal pressure varies with 
temperature at fixed density. The results given below are 
in complete agreement with these expectations. 

We evaluated Eq. ([7]) by performing MD simulations 
of water in the canonical ensemble for various volumes at 
T = 300 K. The dependence of the internal pressure on 
density at T = 300 K is shown in Fig. [71 Note that the 
total internal pressure. Pi, becomes increasingly negative 
as p is increased, in direct opposition to the ap'^ depen- 
dence given by the vdW equation of state. However, it 
is known that as the density of a LJ fluid is increased to 
high values so that neighboring repulsive cores begin to 



overlap, the total Pi exhibits a maximum, after which the 
internal pressure becomes increasingly negative from the 
dominant contribution of the repulsive interactions [l3| . 

As shown in Fig. [2] there is substantial overlap of the 
repulsive LJ cores between nearest neighbors in SPC/E 
water. The repulsive interactions from these LJ cores 
dominate the density dependence of both e and Pi for 
SPC/E and related water models, as evidenced by the 
similarity of the internal pressures of both the full and 
GTRC water models in Fig. H Although e"? > e^'' for 
all density, e* does not exhibit very large changes upon 
increasing density, a direct consequence of the ability of 
water to maintain its hydrogen bond network under the 
conditions studied. Thus the density dependence of the 
internal pressure of SPC/E water is actually similar to 
that of a L J fluid but one at a very high effective density 
with substantial overlap of neighboring cores. 

In addition to the anomalous density dependence of 
Pi , the temperature dependence of the internal pressure 
of water has also been called an anomaly [33 • For most 
organic liquids (and vdW fluids), the internal pressure 
decreases with increasing temperature, but that of water 
increases when the temperature is increased, as shown in 



Fig. [8l Using the concepts presented above, we can ra- 
tionalize this behavior in terms of molecular interactions. 
By decomposing Pi into its electrostatic and LJ compo- 
nents, we find that P' dominates the temperature depen- 
dence of the internal pressure, increasing with increasing 
temperature, while P/"^ is dominated by repulsive in- 
teractions at all temperatures studied, as evidenced by 
its negative value for all T. As the temperature of the 
system is increased, the number of ideally tetrahedrally 
coordinated water molecules decreases, and the hydrogen 
bond network becomes increasingly "flexible" . Therefore, 
if one increases the volume of the system at high T, wa- 
ter will more readily expand to fill that volume. But an 
increase in the electrostatic energy will also occur due to 
a slight decrease in the number of (favorable) hydrogen 
bonding interactions. This will happen to a lesser extent 
at low temperatures, when the hydrogen bond network is 
more rigid and the thermal expansivity of water is lower. 



IV. UNBALANCED FORCES IN 

NONUNIFORM AQUEOUS MEDIA FROM THE 

VIEWPOINT OF LMF THEORY 

In contrast to uniform systems, a net cancellation of 
long ranged forces does not occur in nonuniform envi- 
ronments, and these unbalanced forces can cause signif- 
icant changes in the structure and thermodynamics of 
the system [2^ [S^- As shown above, the bulk struc- 
ture of both the GT and GTRC models are very similar 
to that of the full water model at a given temperature 
and density. But interfacial structure and coexistence 
thermodynamic properties of the uncorrected reference 
models can be very different. For example, GTRC wa- 
ter still has a self-maintained liquid-vapor (LV) interface 
at T = 300 K as illustrated in Fig. [9l even though the 
LJ attractions are ignored, because of the strong charge 
pairing leading to hydrogen bond formation. However 
its 90-10% interfacial width increases to w ~ 4.9 A from 
the w w 3.5 A seen in both GT water and the full water 
model, and the coexisting liquid density of GTRC water 
is about about 15% lower. In contrast, the density pro- 
file of the GT model with LJ interactions fully accounted 
for is in very good qualitative agreement with that of the 
full model. This strongly suggests that if local hydrogen 
bonding is properly taken into account, the equilibrium 
structure of the LV interface of water is governed mainly 
by long ranged LJ attractions, with long ranged dipole- 
dipole interactions playing a much smaller role. It is the 
exact balance of these long ranged interactions we seek 
to examine in this section. 

LMF theory provides a framework in which the aver- 
aged effects of long ranged forces are accounted for by 
an effective external field [SO]- It has previously been 
used mainly as a computational tool to permit very ac- 
curate determination of properties of the full nonuniform 
system while using a numerical simulation of the short 
ranged reference system in the presence of the effective 
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FIG. 9. (Color online) Density profiles of oxygen sites at the 
liquid-vapor interface of Full, GT, and GTRC water models. 
The Gibbs dividing surface of each interface is located at z — 
0. 



field [Tfl [ij, 123 • But the effective or renormalized field 
also gives a convenient and natural measure of the impor- 
tance of long ranged forces in different environments. In 
this section we use the renormalized external fields deter- 
mined directly from simulations of interfaces in the full 
SPC/E water model along with simulations of truncated 
water models to quantitatively examine the relative infiu- 
ence of the local hydrogen bond network and unbalanced 
long-ranged Coulomb and van der Waals forces. 



A. Water-vapor and water-solid interfaces 

We first consider the LV interfaces of the SPC/E, GT, 
and GTRC water models shown in Fig. [HI The removal 
of long-ranged electrostatics in the GT model leaves 
the density distribution virtually unchanged, whereas re- 
moval of the LJ attractions in GTRC water has a sub- 
stantial impact on p{z). To understand this behavior, 
we focus our attention on the impact of the averaged un- 
balanced forces from the long-ranged electrostatic and LJ 
interactions, as determined in LMF theory from the effec- 
tive external fields Vri and cj)]^-^ , respectively and defined 
below. The unbalanced force J-" acting on an oxygen site 
from the LMF potentials is given by 



J-o(r) = -Vr^i^Kr) - goVrVRi(r) 



(8) 



Here go is the partial charge on the oxygen site and 
VRi(r) is the slowly- varying part of the effective elec- 
trostatic field, given by 



VRi(r) 



drV'^(r')ui (|i 



(9) 



where p''{r) is the total charge density of the system. 
The other contribution 0{^'}(r) is the field arising from 



the unbalanced LJ attractions on the oxygen site (where 
the LJ core is centered), given by 



^]^,{r)^ dr'[p{r')-p^]u,{\r~r'\), 



(10) 



with p(r) indicating the nonuniform singlet density distri- 
bution of oxygen sites and p^ defined as the bulk density 
of oxygen sites at the state point of interest [3(|,!39|- Since 
the hydrogen sites lack LJ interactions, the unbalanced 
LMF force acting on a hydrogen site is due exclusively 
to electrostatics, 



J-H(r) = -(?HVrVRi(r). 



(11) 



Given its importance in the density distribution of wa- 
ter, it may seem natural to examine the components of 
the LMF force on the oxygen sites, To{z), shown in the 
inset of Fig. [TUk . Naive examination of the relative mag- 
nitude of these force functions would lead to the con- 
clusion that long-ranged electrostatics are the dominant 
unbalanced force at the LV interface. However, Vri also 
interacts with hydrogen sites and one should instead con- 
sider the net forces from long ranged Coulomb and LJ 
interaction felt by an entire water molecule at these in- 
terfaces. 

This ensemble averaged net molecular force (J-") 
fFig.llOb) clearly indicates that the net unbalanced force 
at an interface is almost entirely due to long-ranged LJ 
attractions from the bulk, which pull water molecules 
away from the interface. The long-ranged Coulomb con- 
tributions to the average force on a water molecule arc es- 
sentially negligible in comparison. This is not surprising 
since water molecules are neutral, and it has previously 
been shown that the small net interfacial electrostatic 
force simply provides a slight torque on water molecules 
in this region [23 ■ This torque has little effect on the 
oxygen density distribution, as illustrated by the good 
agreement of the GT model density profile with that of 
the full water model in Fig. [HI However, it plays a key 
role in determining electrostatic and dielectric properties, 
which are strongly affected by the behavior of the total 
charge density, and here the uncorrected GT model gives 
very poor results [23, [SO] ■ 

It is also instructive to compare the unbalanced long 
ranged forces at the LV interface to those at the liquid- 
solid (LS) interface between water and a model hy- 
drophobic 9-3 LJ wall introduced by Rossky and cowork- 
ers [22], as shown in Fig. [TUb . Despite the large differ- 
ences in the density profiles shown in Fig. llOb . the net 
unbalanced forces on molecules at the LV and LS in- 
terfaces are remarkably similar for all z until molecules 
encounter the harsh repulsion of the wall and an accurate 
sampling of {T{z)) by simulation cannot be made. Wa- 
ter molecules can sample all regions in the liquid-vapor 
interface, leading to a smooth {J-{z)) at smaller z. 

Indeed, the net molecular force due explicitly to a 
configurational average of the attractive ui (r) acting on 
molecules present at each z-position is in outstanding 
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FIG. 10. (Color online) (a) Density distributions as a function 
of the ^-coordinate for the hydrophobic LS interface and the 
LV interface of water, (b) Ensemble averaged net force on a 
water molecule due to Vri (open symbols) and (^^j (closed 
symbols) at the LV (circles) and LS (squares) interfaces and 
the density profiles of the two systems. Solid lines indicate 
the net force due to ui. The black dashed line at z = A 
indicates the position of the hydrophobic wall. The Gibbs 
dividing interface of the LV system is located at z — 2.34 A, 
in order to make comparison with the water-wall interface. 
Inset: Forces on oxygen sites only, determined by evaluating 
the gradient of the corresponding LMF potentials. Labeling 
for the inset is the same as that in (b). 



quantitative agreement with that arising from 0p}(z) for 
all adequately sampled regions in the liquid, as illustrated 
by the solid lines in Fig.llOb. This serves largely as confir- 
mation of the validity of the mean-field treatment inher- 
ent in LMF theory within the liquid slabs. Deviations 
between the two quantities for distances less than the 
Gibbs dividing surface are a reflection of the increasing 
effect of larger force fluctuations due to long- wavelength 
capillary waves not well described by mean field theory. 
The relative magnitudes of the components of (J^) for 
the LV and LS interfaces are strikingly similar, with the 
net unbalanced LJ force (j^ (z; ^j^^J)) reaching its max- 
imum value of slightly less than k-^T / A near the Gibbs 
dividing interface and the repulsive boundary of the wall, 
respectively. 

The similarities of the unbalanced forces at the LV and 
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the hydrophobic LS interfaces of water and the dom- 
inance of the LJ attractions are completely consistent 
with the analogies commonly drawn between these two 
systems [3, lo l ISSl an d used in the LCW theory of hy- 
drophobicity l2J, lla, [S^- A common criticism of LCW 
theory is its apparent neglect of the hydrogen bond net- 
work of water and the use of a van der Waals like expres- 
sion for the unbalanced force at an interface. Although 
some features of the network are implicitly captured by 
using the experimental surface tension and radial dis- 
tribution function of water as input to the theory, elec- 
trostatic effects at the interface, including dipole-dipole 
interactions, are ignored. However, this assumption is 
justified since the averaged effects of long-ranged dipole- 
dipole interactions, accounted for by Vri, are shown to 
indeed be negligible at a hydrophobic interface (Fig. [TU]) . 
LCW theory correctly describes the unbalanced LJ at- 
tractions from the bulk, which dominates the behavior 
at both the liquid-vapor and extended hydrophobic in- 
terfaces. 



B. Crossover from small to large length scale 
hydration 

LCW theory combined the idea of unbalanced forces 
with experimental data in water to predict a crossover in 
the solvation of spherical hydrophobic solutes occurring 
at a radius of about 1 nm [2J| . While the local hydrogen 
bond network can be maintained around smaller solutes, 
some bonds must be broken on larger length scales, lead- 
ing to an enthalpically dominated regime [3, M, ^M ■ Here 
we use LMF theory and a direct analysis of the unbal- 
anced long ranged forces to provide further physical in- 
sight into how this crossover comes about. We study the 
hydration of spherical hydrophobic solutes, which inter- 
act with the oxygen site of water via a solute-water LJ 
potential Usw{r) with a fixed well depth of Ssw = 0.19279 
kcal mol^^ and varying solute- water interaction length- 
scales, CTsw, ranging from 2.0Atol5.0A. This size range 
spans the crossover from small to large length-scale hy- 
drophobic hydration as determined by both simulation 
and theory [3, UM ■ 

As shown earlier, the unbalancing force from LJ attrac- 
tions dominates that due to long ranged electrostatics at 
nonpolar interfaces so we focus only on the LJ forces in 
our analysis of the length scales of hydration. The total 
renormalized external solute-solvent field is the sum of 
the bare solute-solvent field and the slowly-varying LMF 
potential ^^'{(r) in Eq. (fTU)) . which accounts for unbal- 
anced LJ forces from the nonuniform water (oxygen) den- 
sity distribution around the solute. It can be written as 



0i^:'(r)=Msw(r) + 0k-lW 



(12) 



We compare these potentials for small and large LJ so- 
lutes in Fig. [TT] At small solute sizes, the renormalized 
field <j>^ exhibits a repulsive core nearly identical to that 
of the bare Usw, and the effective attractions are hardly 
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FIG. 11. (Color online) Comparison of renormalized solute- 
water vdW potentials for aaw = 2.0 A (squares) and (Ts„ = 
15.0 A (circles). The full solute- water LJ potential Usw is also 
shown for comparison (black dash-dot line). Note that the 
X-axis is scaled by the solute-water interaction length-scale 



altered upon renormalization of the potential, indicat- 
ing only a small unbalanced force around a small solute. 
The drive to maintain the hydrogen bond network around 
small solutes dominates the water structure, and solute- 
water and water-water LJ attractions are found to have 
little effect on the solvation structure Q. 

In contrast, water cannot completely preserve its hy- 
drogen bond network at an interface around a large solute 
and one hydrogen bond per interfacial water molecule 
tends to be broken on average 0, M, [22|, [S^l ■ This cre- 
ates a soft fluctuating interface for which partial drying 
can be induced by unbalanced attractions from the bulk. 
These are re-expressed in LMF theory as an effective re- 
pulsion pushing water away from the solute surface, as 
illustrated by the renormalized potential for the large so- 
lute in Fig. 111! In order for water to wet such an extended 
solute surface, this large effective repulsion must be over- 
come. In hydrophilic surfaces this typically arises from 
strongly attractive polar groups on the surface, such as 
hydrogen bonding sites or charged groups, which also can 
strongly perturb local hydrogen bond configurations. 

To illustrate the interplay between unbalanced forces 
from the bulk, water-water hydrogen-bonding, and the 
length-scale dependence of hydration, we study the hy- 
dration of a small ((Ts„ = 2.0 A) and large (ctsw = 15.0 A) 
LJ solute by SPC/E water and the short-ranged GTRC 
water model. Since the GTRC model lacks LJ attrac- 
tions, the unbalancing force from the bulk is absent but 
the local hydrogen bonding network of water is accurately 
captured. The GTRC model thus properly describes the 
dominant crossover behavior of retaining or breaking hy- 
drogen bonds in the small or large length scale regimes, 
respectively, but the subsequent large length scale inter- 
face properties will be incorrect. Comparing GTRC wa- 
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FIG. 12. (Color online) Singlet density distributions of water oxygen sites, p{r), around LJ solutes with asw = 2.0 A (a) and 
ffsw = 15.0 A (b) for Full and GTRC water models at pressures of P = 1.0 atm and Po = P — Pi, respectively. Data for GTRC 
water in the presence of the renormalized LMF solute potential is also shown for the large solute (GTRC-LMF). 



tcr with the full water model provides more insight into 
the relative importance of the local network and the un- 
balancing force on the structuring of water around each 
solute. 



Singlet density distributions of the water oxygen sites 
around the small and large LJ solute are shown in 
Fig.[T2k and[T2b. respectively. In the small solute regime, 
p{r) for GTRC water is essentially identical to that of the 
full water model. This provides dramatic confirmation 
of the standard physical picture that small length scale 
hydrophobicity is almost completely dominated by the 
need to maintain local hydrogen bonding around the so- 
lute. Unbalanced forces from either the LJ or long-ranged 
Coulomb interaction play essentially no structural role in 
this regime. 

In contrast, although the GTRC model correctly de- 
scribes the necessary breaking of local hydrogen bonds 
around a large solute in Fig. [T2b the details of the result- 
ing interface profile are very different from that of the full 
water model. Removal of L J attractions in GTRC water 
eliminates the large effective repulsion at the extended 
solute surface, and uncorrected GTRC water appears to 
wet the surface of the nonpolar solute. However, the in- 
terface is soft and when the unbalanced force from 0^1 
in Fig. [To] as given by LMF theory is also taken into 
account, depletion and the correct structuring of water 
at the solute surface is very accurately described by the 
GTRC-LMF curve in Fig. [Hj. 



on water due to the renormalized LJ solute external field, 



Fs 



^k'] 



drp{r) 



9Msw(r) 



dr 



drpir)- 



dr 



'S [Us 



Fs 



^ki] 



(13) 

As illustrated in Fig. [T31 the unbalanced force due to 
water-water LJ attractions dominates in the large-scale 
regime, and can only be overcome by unphysically large 
solute-water LJ attractions. In the small scale regime, 
solute- water attractions are comparable in magnitude to 
the unbalancing force, and can be larger for certain so- 
lute sizes (the relative magnitudes are dependent upon 
the value of Esw)- However, as the GTRC model shows, 
the important physics in this regime is simply maintain- 
ing the hydrogen bond network, and this imposes a near 
constant solvation structure as the water-solute attrac- 
tions are varied 0, [l3| • 

Upon normalizing the mean unbalanced force by the 
surface area of the solute. 



fs 



^li] 



Fs [0^i] /^nal 



(14) 



a transition from scaling of Fs with solute volume to scal- 
ing with surface area occurs at the small to large length 
scale solvation crossover length of roughly 10 A, as ev- 
idenced by the plateau in fs depicted in the inset of 
Fig. 1131 This plateau is similar to that which occurs in 
the solvation free energy [13, l2Jl , and provides another 
indicator of the transition from small to large length scale 
hydrophobicity. 



Although the qualitative dependence of the unbalanced 
force on solvation length scale is depicted in Fig. [Til a 
more quantitative metric of this behavior is desirable. To 
that end, we introduce the mean solvation force acting 



V. CONCLUSIONS 

In this work, we have examined the different roles of 
short and long ranged forces in the determination of the 
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FIG. 13. (Color online) Components of the mean solvation 
force Fs due to 0p{ and itsw as a function of the solute-water 
interaction length scale, ctsw The inset displays the solvation 
force due to (/ijjj scaled by the surface area of the solute. 



structure and thermodynamics of uniform and nonuni- 
form aqueous systems, using concepts inherent in classi- 
cal perturbation and LMF theory. In particular, we have 
evaluated individually for SPC/E water the contributions 
of (i) all the strong short ranged repulsive and attrac- 
tive interactions that lead to the local hydrogen-bonding 
network, (ii) longer ranged dispersive LJ attractions be- 
tween molecules, and (iii) long ranged dipole-dipole in- 
teractions, and demonstrated a hierarchical ordering of 
their importance in determining several properties of wa- 
ter in uniform and nonuniform systems. 

All of our truncated models accurately describe the lo- 
cal hydrogen bonding network, and as expected, this net- 
work alone is sufficient to match bulk structure as well as 
solvent structure around small hydrophobic solutes pro- 
vided that the bulk density and temperature are accu- 
rately prescribed. Furthermore, the anomalous tempera- 
ture and density dependence of the "internal pressure" of 
water is found to be dominated by the competing short- 
ranged repulsive and attractive forces determining the 
local hydrogen bonding network as well. 

But local network concepts alone cannot capture 
all the complexities of even the simple SPC/E water 
model. While the dispersive LJ attractions between wa- 
ter molecules primarily provide a uniform cohesive energy 
in bulk systems, they strongly influence the structure 
and density proflle of large scale hydrophobic interfaces. 
Their importance provides further support for analogies 
between water at extended hydrophobic interfaces and 
the liquid-vapor interface, and the unbalanced LJ force 
can be used to quantify the transition between small and 
large scale hydrophobicity for simple solutes. 

Although the long-ranged dipolar interactions between 
molecules have only small effects on most of the interfa- 
cial density properties considered here, we have shown 



elsewhere that they are crucial in determining dielectric 
properties of both bulk and nonuniform water. Indeed, 
as will be discussed elsewhere, we have found that elec- 
trostatic quantities may in fact be a sensitive structural 
probe of hydrophobicity in general environments |27| . 

This interaction hierarchy, wherein strong short- 
ranged local interactions alone determine structure in 
uniform environments while the longer ranged forces are 
needed as well to capture other properties could prove 
quite useful in rcflning simple site-site water models. 
Current water models incorporate a vast amount of clever 
engineering and empirical fine-tuning and manage to re- 
produce a variety of different properties through a com- 
plex balance of competing interactions with simple func- 
tional forms. Changes in the potential that improve one 
property generally speaking produce poorer results for 
several others. 

One promising route to a more systematic procedure 
may be sensitivity analysis, in which small perturbations 
of potential parameters are made and the correlated re- 
sponse of a variety of physical observables is quantified. 
By perturbing the relative magnitudes of short and long 
ranged interactions, lordanov et al. found that thermody- 
namic properties of bulk water are most sensitive to small 
changes in the L J repulsions and the short ranged electro- 
static interactions [19|, in agreement with our findings. 
A new water model was then proposed by optimizing 
parameters to reproduce a specific bulk thermodynamic 
quantity (the internal energy) in an attempt to correct 
the deficiencies present in a previously developed water 
potential. 

However, the theoretical scheme of splitting the poten- 
tial described in this paper may provide a more concrete 
and physically suggestive path to incrementally match 
various known physical quantities for water without ru- 
ining the fitting of previous quantities, and one could 
combine an approach like sensitivity analysis with the 
conceptual framework presented herein to systematically 
optimize a specific water model. 

In particular, it has recently been suggested that the 
accuracy with which a water model can predict the exper- 
imental Tmd correlates well with the accuracy that the 
same model displays in predicting the thermodynamics 
of small-scale hydrophobic hydration Q • Arguably, the 
least justified feature of current simple water models like 
SPC/E is the functional form of the core LJ potential 
uo{r), especially at the very short separations relevant 
for describing local hydrogen bonding as illustrated in 
Fig. [21 One could try to fine-tune a GTRC-type model 
through alteration of the local hydrogen bond network by 
varying the form of the repulsive core in order to match 
the experimental density maximum, as well as other bulk 
properties like the internal pressure, in order to obtain 
a short-ranged system that yields accurate bulk prop- 
erties. Although a detailed discussion of this process is 
beyond the scope of this article, one could try to use some 
type of o ptim ization procedure to determine such poten- 
tials [1^ [2y, I43 ■ Perhaps first principles DFT Simula- 
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System 



iV 



equil ''run 



Bulk Water 
Bulk Water 
LV Interface 
LS Interface 
Small LJ Solutes (ctsw 
Large LJ Solutes (ctsw 



256 220-340 K 5 ns 5 ns 

1000 220-300 K 2 ns 2 ns 

1728 298 K 1 ns 500 ps 

2468 298 K 1 ns 500 ps 

< 10 A) 1000 298 K 2 ns 2 ns 

> 10 A) 6000 298 K 1 ns 1 ns 



TABLE I. Details of the MD simulations performed in this 
work. A'' and T refer to the number of water molecules and 
the temperature of the system, respectively. The equilibration 
times and data collection times are denoted by teq^u and trun, 
respectively. 



tions [211, [33 could be used to provide a more fundamen- 
tal description of the local netvifork. Subsequently, the 
structure and thermodynamics of nonuniform systems, 
which require dispersions and long ranged Coulomb in- 
teractions, could be used to parametrize the long-ranged 
interactions. 
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Appendix: Simulation Details 

All molecular dynamics simulations were performed us- 
ing modified versions of the DL_POLY software pack- 
age 3J| and the SPC/E water model Q or its variants 
described in Section 2. The equations of motion were 
integrated using the leapfrog algorithm with a timestep 
of 1 fs m while maintaining constant temperature and 
pressure conditions through the use of a Berendsen ther- 
mostat and barostat respectively Q. 



1. Bulk water simulations 

The evaluation of electrostatic interactions in bulk sim- 
ulations of the full SPC/E water model employed the 
standard Ewald summation method using a real space 
cutoff of 9.5 A, unless this was larger than half of the box 
length, in which case the cutoff was set to half of the box 



length [l|. Short-ranged electrostatic interactions in the 
GT and GTRC reference systems, as well as LJ interac- 
tions in all systems, were truncated at the real space cut- 
off length used in the analogous full system. Simulations 
of bulk water were performed with N ~ 1000 molecules 
in the isothcrmal-isobaric (NPT) ensemble to determine 
the density maximum and with N = 256 molecules in 
the canonical (NVT) ensemble to determine P{T) and 
the internal pressure. The duration of equilibration and 
production runs, as well as the temperatures sampled are 
listed in Table 1. The internal pressure in Eq. ([7]) was 
calculated by evaluating s{v) for numerous values of v 
at each T. The function e{v) was then fit to a polyno- 
mial, which was differentiated at the desired v to yield 
the internal pressure. 



2. Simulation of nonuniform systems 

In order to generate starting configurations for the LV 
and LS interfacial systems discussed in Section 5, we 
first equilibrated N water molecules in a cubic geome- 
try, where N is listed in Table 1. The z-dimcnsion of the 
system was then elongated to more than three times the 
X- and y-dimensions, and in the case of the LS interface, 
a wall potential of the form 



U^,{z) 



B 



I ^ ^W I I ^ ^W I 



(A.l) 



was added at Zj^, = and the parameters A and B are 
given in Ref. [23|. In order to ensure water molecules 
did not approach the wall from z < 0, a repulsive wall 
was added at large z to constrain the water molecules 
to the desired region of the simulation cell while still al- 
lowing a large vacuum region for the formation of a va- 
por phase. Electrostatic interactions were handled using 
the corrected Ewald summation method for slab geome- 
tries [4^ with a real space cutoff of 11.0 A. which was also 
the cutoff distance for LJ and short-ranged electrostatic 
interactions. 

Molecular dynamics simulations of the hydration of LJ 
solutes were performed in the NPT ensemble using stan- 
dard Ewald summation to evaluate the electrostatic in- 
teractions, with a real space cutoff of 11.0 A. The short 
ranged vq potential and the water-water LJ potential 
were also truncated at 11.0 A. The LJ solute was rep- 
resented by a fixed external potential centered at the ori- 
gin, and the solute-water interactions were truncated at 
one-half the length of the simulation cell. The number of 
molecules and simulation times are listed in Table 1. 
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